Celem projektu jest określenie jakie mogą być główne przyczyny stopniowego zmniejszania się długości śledzi oceanicznych wyławianych w Europie.
Zbiór danych zostanie wczytany z pliku CSV, następnie musi zostać poddany wstępnemu oczyszczaniu. # Opis zbioru danych Analiza dotyczy zbióru danych na temat połowu śledzia oceanicznego w Europie. Do analizy zebrano pomiary śledzi i warunków w jakich żyją z ostatnich 60 lat. Dane były pobierane z połowów komercyjnych jednostek. W ramach połowu jednej jednostki losowo wybierano od 50 do 100 sztuk trzyletnich śledzi.
Poniżej znajdują się szczegółowe opisy konkretnych atrybutów:
| Nazwa kolumny | Opis | Dodatkowa Informacja |
|---|---|---|
| length | długość złowionego śledzia | [cm] |
| cfin1 | dostępność planktonu | [zagęszczenie Calanus finmarchicus gat. 1] |
| cfin2 | dostępność planktonu | [zagęszczenie Calanus finmarchicus gat. 2] |
| chel1 | dostępność planktonu | [zagęszczenie Calanus helgolandicus gat. 1] |
| chel2 | dostępność planktonu | [zagęszczenie Calanus helgolandicus gat. 2] |
| lcop1 | dostępność planktonu | [zagęszczenie widłonogów gat. 1] |
| lcop2 | dostępność planktonu | [zagęszczenie widłonogów gat. 2] |
| fbar | natężenie połowów w regionie | [ułamek pozostawionego narybku] |
| recr | roczny narybek | [liczba śledzi] |
| cumf | łączne roczne natężenie połowów w regionie | [ułamek pozostawionego narybku] |
| totaln | łączna liczba ryb złowionych w ramach połowu | [liczba śledzi] |
| sst | temperatura przy powierzchni wody | [°C] |
| sal | poziom zasolenia wody | [Knudsen ppt] |
| xmonth | miesiąc połowu | [numer miesiąca] |
| nao | oscylacja północnoatlantycka | [mb] |
library(knitr)
library(ggplot2)
library(polycor)
library(heatmaply)
library(tidyr)
library(plotly)
library(VIM)
library(caret)
library(klaR)
library(dplyr)
set.seed(23)
raw_data <- read.csv(file= "sledzie.csv", header= TRUE, sep= ",", na.strings= "?")
str(raw_data)
## 'data.frame': 52582 obs. of 16 variables:
## $ X : int 0 1 2 3 4 5 6 7 8 9 ...
## $ length: num 23 22.5 25 25.5 24 22 24 23.5 22.5 22.5 ...
## $ cfin1 : num 0.0278 0.0278 0.0278 0.0278 0.0278 ...
## $ cfin2 : num 0.278 0.278 0.278 0.278 0.278 ...
## $ chel1 : num 2.47 2.47 2.47 2.47 2.47 ...
## $ chel2 : num NA 21.4 21.4 21.4 21.4 ...
## $ lcop1 : num 2.55 2.55 2.55 2.55 2.55 ...
## $ lcop2 : num 26.4 26.4 26.4 26.4 26.4 ...
## $ fbar : num 0.356 0.356 0.356 0.356 0.356 0.356 0.356 0.356 0.356 0.356 ...
## $ recr : int 482831 482831 482831 482831 482831 482831 482831 482831 482831 482831 ...
## $ cumf : num 0.306 0.306 0.306 0.306 0.306 ...
## $ totaln: num 267381 267381 267381 267381 267381 ...
## $ sst : num 14.3 14.3 14.3 14.3 14.3 ...
## $ sal : num 35.5 35.5 35.5 35.5 35.5 ...
## $ xmonth: int 7 7 7 7 7 7 7 7 7 7 ...
## $ nao : num 2.8 2.8 2.8 2.8 2.8 2.8 2.8 2.8 2.8 2.8 ...
Zbiór zawiera 52582 rekordów rozmieszczonych w 16 kolumnach (z czego jedna jest kolumną porządkową).
Zmienna xmonth, która reprezentuje miesiąc połowu powinna zostać zamieniona z ciągłej na kategoryczną, by nie traktować jej jako liczbę. Zmienna totaln, która reprezentuje łączną liczbę ryb złowionych, powinna zostać zmieniona na całkowitą.
raw_data <- raw_data %>%
mutate(xmonth= as.factor(xmonth), totaln= round(totaln), totaln= as.integer(totaln))
#Puste dane liczbowo na kolumnę
apply(raw_data, 2, function(x){ sum(is.na(x)) })
## X length cfin1 cfin2 chel1 chel2 lcop1 lcop2 fbar recr
## 0 0 1581 1536 1555 1556 1653 1591 0 0
## cumf totaln sst sal xmonth nao
## 0 0 1584 0 0 0
#Puste dane procentowo na kolumnę
apply(raw_data, 2, function(x){ sum(is.na(x)) / length(x) })
## X length cfin1 cfin2 chel1 chel2
## 0.00000000 0.00000000 0.03006732 0.02921152 0.02957286 0.02959188
## lcop1 lcop2 fbar recr cumf totaln
## 0.03143661 0.03025750 0.00000000 0.00000000 0.00000000 0.00000000
## sst sal xmonth nao
## 0.03012438 0.00000000 0.00000000 0.00000000
Zbiór zawiera również wartości puste - te pojawiają się głównie w kolumnach z informacją o dostępności planktonu oraz temperaturze przy powierzchni wody.
aggr(raw_data, plot= TRUE,
col= c('#fa9fb5', '#2b8cbe'),
numbers= TRUE,
prop= FALSE,
bars= FALSE,
labels= names(raw_data),
cex.axis= 0.8,
ylab=c("Histogram brakujących danych","Wzorzec"))
Jak zobrazowano na wykresie powyżej, rozkład wartości pustych w kolumnach:
Usuwając wiersze z wartością NA, utracilibyśby stosunkowo dużo danych - lepszym pomysłem jest zastąpienie wartości brakującej średnią z konkretnego połowu. Bazując na fakcie, iż kolumny totaln, xmonth, nao definiują konkretny połów oraz nie zawierają one żadnych wartości pustych, posłużą one do grupowania. Dane zostały zgrupowane względem połowów, a następnie wartości puste zostały zamienione na średnią z tych połowów.
data <- raw_data %>%
group_by(totaln, xmonth, nao) %>%
mutate_each(funs(replace(., which(is.na(.)),
mean(., na.rm=TRUE))))
no_x <- data %>% select(-X)
sum(duplicated(no_x))
## [1] 45694
Możemy zaobserwować, iż 45694 rekordów to duplikaty. Pojawiają się one wewnątrz jednego połowu, dlatego usunięcie ich nie wpłynie negatywnie, ani nie sfałszuje danych. Dla uproszczenia grafów oraz dalszych obliczeń, wszystkie duplikaty zostały usunięte, tym samym zbiór danych uszczuplił się do 6888 rekordów.
w_duplicates <- unique(no_x[, 1:15])
w_duplicates <- w_duplicates %>%
ungroup %>%
mutate(X = seq_len(n())) %>%
select(X, everything())
str(w_duplicates)
## Classes 'tbl_df', 'tbl' and 'data.frame': 6888 obs. of 16 variables:
## $ X : int 1 2 3 4 5 6 7 8 9 10 ...
## $ length: num 23 22.5 25 25.5 24 22 23.5 22.5 22 24.5 ...
## $ cfin1 : num 0.0278 0.0278 0.0278 0.0278 0.0278 ...
## $ cfin2 : num 0.278 0.278 0.278 0.278 0.278 ...
## $ chel1 : num 2.47 2.47 2.47 2.47 2.47 ...
## $ chel2 : num 21.4 21.4 21.4 21.4 21.4 ...
## $ lcop1 : num 2.55 2.55 2.55 2.55 2.55 ...
## $ lcop2 : num 26.4 26.4 26.4 26.4 26.4 ...
## $ fbar : num 0.356 0.356 0.356 0.356 0.356 0.356 0.356 0.356 0.356 0.356 ...
## $ recr : num 482831 482831 482831 482831 482831 ...
## $ cumf : num 0.306 0.306 0.306 0.306 0.306 ...
## $ totaln: int 267381 267381 267381 267381 267381 267381 267381 267381 267381 267381 ...
## $ sst : num 14.3 14.3 14.3 14.3 14.3 ...
## $ sal : num 35.5 35.5 35.5 35.5 35.5 ...
## $ xmonth: Factor w/ 12 levels "1","2","3","4",..: 7 7 7 7 7 7 7 6 6 6 ...
## $ nao : num 2.8 2.8 2.8 2.8 2.8 2.8 2.8 2.8 2.8 2.8 ...
Zbiór danych po oczyszczaniu zmniejszył się do 6888 wierszy, liczba kolumn pozostała niezmieniona i wynosi 16. ## Analiza arybutów
knitr::kable(summary(w_duplicates))
| X | length | cfin1 | cfin2 | chel1 | chel2 | lcop1 | lcop2 | fbar | recr | cumf | totaln | sst | sal | xmonth | nao | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| Min. : 1 | Min. :19.00 | Min. : 0.00000 | Min. : 0.0000 | Min. : 0.000 | Min. : 5.238 | Min. : 0.3074 | Min. : 7.849 | Min. :0.0680 | Min. : 140515 | Min. :0.06833 | Min. : 144137 | Min. :12.77 | Min. :35.40 | 8 : 879 | Min. :-4.89000 | |
| 1st Qu.:1723 | 1st Qu.:24.00 | 1st Qu.: 0.02778 | 1st Qu.: 0.2500 | 1st Qu.: 2.469 | 1st Qu.:15.030 | 1st Qu.: 2.5479 | 1st Qu.:20.094 | 1st Qu.:0.1580 | 1st Qu.: 364794 | 1st Qu.:0.11008 | 1st Qu.: 307276 | 1st Qu.:13.64 | 1st Qu.:35.51 | 10 : 879 | 1st Qu.:-1.69000 | |
| Median :3444 | Median :25.50 | Median : 0.14158 | Median : 0.3714 | Median : 4.811 | Median :21.435 | Median : 5.9167 | Median :24.859 | Median :0.3320 | Median : 459347 | Median :0.21476 | Median : 539558 | Median :13.98 | Median :35.51 | 7 : 747 | Median : 0.20000 | |
| Mean :3444 | Mean :25.32 | Mean : 0.55913 | Mean : 1.7403 | Mean : 8.801 | Mean :21.157 | Mean : 11.3557 | Mean :27.683 | Mean :0.3202 | Mean : 543028 | Mean :0.21417 | Mean : 523418 | Mean :13.94 | Mean :35.52 | 9 : 680 | Mean : 0.08938 | |
| 3rd Qu.:5166 | 3rd Qu.:26.50 | 3rd Qu.: 0.36032 | 3rd Qu.: 1.5701 | 3rd Qu.: 9.667 | 3rd Qu.:26.324 | 3rd Qu.: 12.4959 | 3rd Qu.:35.153 | 3rd Qu.:0.4250 | 3rd Qu.: 774993 | 3rd Qu.:0.28116 | 3rd Qu.: 763083 | 3rd Qu.:14.21 | 3rd Qu.:35.52 | 6 : 559 | 3rd Qu.: 1.80000 | |
| Max. :6888 | Max. :32.50 | Max. :37.66667 | Max. :19.3958 | Max. :75.000 | Max. :57.706 | Max. :115.5833 | Max. :68.736 | Max. :0.8490 | Max. :1565890 | Max. :0.39801 | Max. :1015595 | Max. :14.73 | Max. :35.61 | 5 : 521 | Max. : 5.08000 | |
| NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | (Other):2623 | NA |
data_dist <- w_duplicates %>%
select(-X) %>%
melt
ggplot(data_dist, aes(x= value)) +
geom_density(fill= "#2b8cbe") +
facet_wrap(~variable, scales= "free") +
theme(axis.text.x = element_text(angle = 90, hjust = 1))
Zmienne, poza length, nie mają rozkładu normalnego.
heatmaply(hetcor(as.data.frame(w_duplicates)), k_col = 2, k_row = 3)
Z powodu różnych klas kolumn, np. xmonth jest zmienną kategoryczną, length ciągłą a X porządkową. Została wyliczona heterogeniczna macierz korelacji, metodą hetcor z biblioteki ploycor.
p <- ggplot(w_duplicates, aes(x= X, y= length)) +
geom_line(alpha= 0.3) +
geom_smooth(method= "gam", formula= y ~ s(x, k= 100), size= 1) +
ggtitle("Zmiana rozmiaru złowionego śledzia w czasie")
ggplotly(p)
Jako, że dane zostały uporządkowane chronologicznie, długość śledzia w czasie prezentowany jest przy użyciu liczby porządkowej X. Z powodu ilości danych, który znacznie obniża czytelność wykresu, została zastosowana metoda smooth, która pozwoli na odkrycie ogólnego wzorca. Użycie wygładzenia liniowego, nie byłoby dostatecznie odpowiednie dla zebranego zestawu danych, dlatego został użyty uogólniony model addytywny gam.
Największa korelacja dotyczy par opisujących planktony, lcop1 i chel1 oraz lcop2 i chel2. Duży współczynnik korelacji pomiędzy cumf oraz totaln, przez co możemy wnioskować, iż wraz ze wzrostem łącznej liczby ryb złowionych w połowie rośnie natężenie połowów. Co więcej możemy zaobreswować korelację pomiędzy cumf oraz fbar - łączne roczne natężenie połowów było wysokie tak samo jak ich intensywność.
Regresor ma za zadanie przewidzieć rozmiary śledzia w kolejnych połowach. Dane zostały podzielone na dwa zbiory: uczący i testowy, z czego 75% całego zbioru zostało potraktowane jako uczące. Uczenie odbyło się przy użyciu metody Repeated Cross-Validation, z powodu niewielkich różnic wartości w zbiorze zastosowano liczbę powtórzeń na poziomie 5 z liczbą powtórzen 2. Model jest tworzony w opariu o model klasyfikacyjny Random Forrest.
fit <- lm(length ~ ., data = no_x)
# Miara R^2
summary(fit)$r.squared
## [1] 0.3315561
# Błąd średnio-kwadratowy
rmse <- function(num) sqrt(sum(num^2)/length(num))
rmse(fit$residuals)
## [1] 1.351338
inTraining <-
createDataPartition(
y = no_x$length,
p = 0.75,
list = FALSE)
training <- no_x[ inTraining,]
testing <- no_x[-inTraining,]
ctrl <- trainControl(
method = "repeatedcv",
number = 2,
repeats = 5)
fit <- train(length ~ .,
data = training,
method = "rf",
trControl = ctrl,
ntree = 2)
## Loading required package: randomForest
## randomForest 4.6-12
## Type rfNews() to see new features/changes/bug fixes.
##
## Attaching package: 'randomForest'
## The following object is masked from 'package:dplyr':
##
## combine
## The following object is masked from 'package:ggplot2':
##
## margin
fit
## Random Forest
##
## 39438 samples
## 14 predictor
##
## No pre-processing
## Resampling: Cross-Validated (2 fold, repeated 5 times)
## Summary of sample sizes: 19719, 19719, 19720, 19718, 19718, 19720, ...
## Resampling results across tuning parameters:
##
## mtry RMSE Rsquared
## 2 1.179349 0.4903080
## 13 1.155392 0.5115404
## 24 1.157489 0.5099190
##
## RMSE was used to select the optimal model using the smallest value.
## The final value used for the model was mtry = 13.
plot(fit)
rfClasses <- predict(fit, newdata = testing)
summary(rfClasses)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 22.09 24.56 25.35 25.31 26.24 28.57
df<-data.frame(rfClasses)
ggplot(df, aes_string(x = rfClasses)) +
geom_histogram(bins= 100, fill= "#0087BD") +
ggtitle("Przewidywany rozmiar śledzia") +
theme_bw() +
labs(x= "Rozmiar śledzi", y="Liczba")
fit_rf <- randomForest(length ~ ., no_x)
fit_rf
##
## Call:
## randomForest(formula = length ~ ., data = no_x)
## Type of random forest: regression
## Number of trees: 500
## No. of variables tried at each split: 4
##
## Mean of squared residuals: 1.319012
## % Var explained: 51.72
importance_df <- importance(fit_rf)
importance_df <- data.frame(var = rownames(importance_df), importance = importance_df[, 1])
importance_df$var <- factor(importance_df$var, levels = importance_df[order(importance_df$importance), "var"])
ggplot(importance_df, aes(x = var, y = importance)) +
geom_bar(stat = "identity", fill = "#2b8cbe") +
ggtitle("Ważność zmiennych") +
theme_bw()
ggplot(data, aes(x = length, y = sst)) +
geom_smooth() +
ggtitle("Zależność długości śledzia od temperatury przy powierzchni wody") +
theme_bw()
## `geom_smooth()` using method = 'gam'
Jak możemy zaobserwować na powyższym wykresie, długość śledzia maleje wraz ze wzrostem temperatury przy powierzchni wody.
ggplot(data, aes(X, sst)) +
geom_smooth() +
theme_bw() +
ggtitle("Zmiana temperatury wody w czasie") +
labs(x= "Czas - l.porzadkowa", y="Temperatura[°C]")
## `geom_smooth()` using method = 'gam'
Natomiast temperatura rosła przez ostatnie lata, co spowodowało znaczne obniżenie długości wyławianych śledzi.
ggplot(data, aes(X, nao)) +
geom_smooth() +
theme_bw() +
ggtitle("Zmiana Oscylacji Północnoatlantyckiej w czasie") +
labs(x= "Czas - l.porzadkowa", y="Oscylacja Północnoatlantycka")
## `geom_smooth()` using method = 'gam'
Konkludując powyższe informacje, możemy postawić diagnozę problemu - w ostatnich latach znacznie wzrosła temperatura przy powierzchni wody, co negatywnie wpłynęło na długość wyławianego śledzia. Wpływ na to ma również zmiana oscylacji północnoatlantyckiej - jest to zjawisko związane z globalną cyrkulacją powietrza i wody oceanicznej, ujawnia się poprzez fluktuacje takich parametrów, jak ciśnienie, temperatura, prędkość wiatru, ilość opadów.